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Abstract. Image reconstruction in X ray tomography consists in determining an 
object from its projections. In many applications such as non destructive testing, we look 
for an image who has a constant value inside a region (default) and another constant 
value outside that region (homogeneous region surrounding the default). The image 
reconstruction problem becomes then the determination of the shape of that region. In 
this work we model the object (the default region) as a polygonal disc and propose a new 
method for the estimation of the coordinates of its vertices directly from a very limited 
number of its projections. 
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1. Introduction 

Gammagraphy is a well known technique in non destructive testing (NDT) and non 
destructive evaluation (NDE) apphcations. Tomographic image reconstruction in 
these apphcations is more recent and consists of determining an object from its 
projections. The relation between the object f{x,y) and its projections p(r, (/>) is 
frequently modeled by the Radon transform: 

p{r,(f>) = jj f {x,y)S{r — X cos (f> — y sin (f)) dx dy (1) 

In many image reconstruction applications, especially in NDT and NDE, we look 
for an image f{x,y) who has a constant value ci inside a region (default region 
D) and another constant value C2 outside that region (homogeneous surrounding 
safe region), e.g. metal & air. The image reconstruction problem becomes then 
the determination of the shape of the default region. 

In this communication, without loss of generality, we assume that ci = 1 and 
C2 = 0: 

" { elsewhere ' ' 

where D represents the default region. 

There has been many works in image reconstruction and computed tomography 
dealing with this problem. To emphasis the originality and the place of this work, 
we give here a summary of the different approaches for this problem: 
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• In the first approach, one starts by discretizing the equation (|l]) to obtain: 

p^Hf + n (3) 

where, / is the discretized values of the object f(x,y) (the pixel values of the 
image), p is values of the projection data p{r,(j)), n is a vector to represent the 
modeling and measurement errors (noise) and H the discretized Radon operator. 
Then the solution is defined as the argument which minimizes the regularization 
criterion 

J{f)^\\p^Hf\\' + Xn{f), (4) 
where A is the regularization parameter. 

i7(/) has to be chosen appropriately to reflect the fact that / must represent a bi- 
nary image. This is the classical approach of general image reconstruction problem. 
In fact, one can also interpret J(/) as the maximum a posteriori (MAP) criterion 
in the Bayesian estimation framework where Q{f) = ||p — HfW^ represents the 
likelihood term and exp [^XQ{f)] the prior probability law. 

This approach has been used with success in many applications (e.g. [|], ||, ^, ^) 
but the cost of its calculation is huge due to the great dimension of /. Many works 
have been done on choosing appropriate regularization functionals or equivalently 
appropriate prior probability laws for / to enforce some special properties of the 
image such as smoothness, positivity or piecewise smoothness [|[ H, ||, ||, [l0|| . 
Among these, one can mention mainly two types of functions for r2(/): 
Entropic laws: 

N 

17(/) = ^0(/j), with 0(a:) = -a;logx,logx, • ■ •} 

Homogeneous Markovian laws: 

^ fx] 

X! ^ifi^fi)^ '^i*^ (l>ix,y) ^ Ux- y^, -\x - 2/1 log-, log cosh I x -y\,---> 

3=1 isAAj y J 

See for example |^ for the entropic laws, 0, ||] for scale invariant markovian laws 
and ^ for other specific choices. 

• In the second approach, one starts by giving a parametric model for the object 
and then tries to estimate these parameters using least squares (LS) or maximum 
likelihood (ML) methods. In general, in this approach one chooses a parametric 
model such as superposition of circular or elliptical discs to be able to relate 
analytically the projections to these parameters. For example, for a superposition 
of elliptical discs we have: 

K 

f{x,y) = Ydkfkix,y) (5) 

k=l 
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with 

f(ry)-i^ if i^-akr + {y-Pk?<gli0), 

JkKx,y} - I Q elsewhere ' 

where 



[6) ^^al cos2 e + bl sin^ 0. (7) 

and where 9 ~ {d^, ak, Pk,ak,bk,k = 1, • • • , K} is a vector of parameters defining 
the parametric model of the image (density values, coordinates of the centers 
and the two diameters of the ellipses). It is then easy to calculate analytically 
the projections and the relation between the data and the unknown parameters 
becomes: 

p{r, 4,) = h{r, (j); 6) + n(r, 4,) (8) 

where 

K 



h{r, ^; 0) = ^ dkPk{r, (f>) with pk[r, (j)) = 



fc=i 



elsewhere 

(9) 

The LS or the ML estimate when the noise is assumed to be Gaussian is then 
given by: 

e = argmin{||p(r,0) - /i(r,0;6')f } (10) 




This approach has also been used with success in image reconstruction [gT| |l^, |13|, 
|l|, |l|. But, the range of applicability of these methods is limited to the cases 
where the parametric models are actually appropriate. 

• In the third approach which is more appropriate to our problem of shape recon- 
struction, one starts by modeling directly the contour of the object by a function, 
say g{9) such as: 

D = {ix,y):p'i9)=x'+y' <g'ie)}. (11) 

The next step is then to relate the projections p(r, 0) to g{9) which, in this case 
is: 

p{r,cj))^ i / 5[r ~ pcos{(i)~e))pApAe. (12) 



and finally to discretize this relation to obtain: 

p^h{g)+n (13) 

where g represents the discretized values of g{9) defining the contour of the object 
and h{g) represents the discretized version of the nonlinear operator (^2|) relating 
projection data p and g. Then, one defines the solution as the argument which 
minimizes 

J{g)^\\p-h{gW+Xn{g), (14) 
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where has to be chosen appropriately to reflect some regularity property of 
the object's contour. 

In this case also one can consider J{g) as the MAP criterion with Q{g) = 
\\p — /i(gr)|p as the likelihood term and ri(g) as the prior one. 
This approach has been used in image restoration , but it seems to be new in 
image reconstruction applications and the proposed method in this work is in this 
category. The originality of our work is to model the contour of the object by a 
piecewise linear function which means that the object is modeled as a polygonal 
disc whose vertices are estimated directly from the projection data. 

Now before going further in details, let compare this last approach with the 
first one, by noting the following: 

• In (H) and (^) , / represents the pixel values of the image (a very great dimensional 
vector depending on the image dimensions), but in ( p^ ) and (p^, g represents the 
discretized values of g{9) defining the contour of the object. The dimension of this 
vector is moderate and independent of the image dimensions. 

In (|) and Hf is a linear function of / and so Q{f) is a quadratic function 
of it, but in ( Jl^ ) and (14), h{g) is not a linear function of g and so Q{g) will not 
be a quadratic function of it. 

We will discuss more the consequences of these remarks in the next section. 



2. Proposed method 

In this paper we propose to model the contour of the object (default region) as 
a periodic piecewise linear function or equivalcntly to model the shape of the 
object as a polygonal disc with a great number N of vertices to be able to 
approximate any shape. Then we propose to estimate directly the coordinates 
{{xj, yj),j = 1, ■ • ■ , N] of the vertices of this polygonal disc from the projection 
data (see Fig. [2]). 

The idea of modeling the shape of the object as a polygonal disc is not new and 
some works have been done in image reconstruction applications, but, in general 
in these works, a hypothesis of convexity of the polygonal disc has been used which 
is very restrictive in real applications. In our work we do not make this hypothesis 
and also we choose N appropriately great to to be able to approximate any shape. 

As we deal with inverse problems, the solution is then defined as the argument 
which minimizes the following criterion 

J{z)^\\p-h{z)\\'' + mz), (15) 

where z = x + iy is a complex vector whose real and imaginary parts represent 
the X and the y coordinates of the polygon vertices, h{z) represents the direct 
operator which calculates the projections for any given z and n{z) is chosen to 
be a function which refiects the regularity of the object contour. In this work we 
used the following: 

N 

1^(2) = ^ - 2zj- + (16) 
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Figure 1: Proposed shape reconstruction modeling. 

Note that \zj-i—2zj+Zj+i p is just the Euclidian distance between the point Zj and 
the line segment passing through zj-i and Zj+i and so this choice favors a shape 
whose local curvature is limited. We can also give a probabilistic interpretation 
to this choice. In fact we can consider Zj as random variables with the following 
Markovian law: 



p{zj\z) = p{zj\zj^i,Zj+i) cx exp 



2z, 



Zj + l\ 



(17) 



Other functions arc possible and are studied in this work. 

In both cases, the criterion J{z) is multimodal essentially due to the fact that 
h(z) is a nonlinear function of z. Calculating the optimal solution corresponding 
to the global minimum of ( |l5| ) needs then carefully designed algorithms. For this 
we propose the following strategies: 

• The first is to use a global optimization technique such as simulated anneal- 
ing (SA). This technique has given satisfactory result as it can be seen from the 
simulations in the next section. However, this algorithm needs a great number of 
iterations and some skills for choosing the first temperature and cooling schedule, 
but the overall calculations is not very important due to the fact that we do not 
need to calculate the gradient of the criterion (|l5|) . 

• The second is to find an initial solution in the attractive region of the global 
optimum and to use a local descent type algorithm to find the solution. 

The main problem here is how to find this initial solution. For this, we used 
a moment based method proposed by Milanfar, Karl & Wilsky Q which 
is accurate enough to obtain an initial solution which is not very far from the 
optimum. The basic idea of this method is to relate the moments of the projections 
to the moments of a class of polygonal discs obtained by an affine transformation 
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of a centered regular polygonal disc, and so to estimate a polygonal disc whose 
vertices are on an ellipse and whose moments up to the second order matches those 
of the projections. 

However, there is no theoretical proof that this initial solution will be in the 
attractive region of the global optimum. In the simulation results section we will 
show some results comparing the performances of these two methods as well as a 
comparison with some other classical methods. 

3. Simulation results 

To measure the performances of the proposed method and keeping the objective 
of using this method for NDT applications where the number of projections are 
very limited, we simulated a case where the object is a polygonal disc with = 40 
vertices (hand-made) and calculated its projections for only 5 directions: 

(j) = {-45, -22.5, 0, +22.5, +45 degrees} 



150- 




Figure 2: Original image and noisy projections. 

Then, we added some noise (white, Gaussian and centered) on them to simulate 
the measurement errors. The S/N ratio was chosen 20dB. Finally, from these data 
we estimated the solution by either of the two proposed methods. Figures |3] and 
^ show these results. 



A BAYESIAN APPROACH TO SHAPE RECONSTRUCTION 



7 



In Fig. |3], we give the reconstruction results obtained by simulated annealing 
(SA) algorithm and in Fig. |3] those obtained by a moment-based initialization and 
a local descent-based optimization algorithm. Note that, the SA is independent 
of initialization, however, in these figures we show the results obtained by the 
proposed method. 



o: orig, +: init, *: rec. ^ ^q* critere J=J1+J2 




-60 -40 -20 20 40 60 10 20 30 40 50 60 70 80 90 100 



X iterations 

Figure 3: Reconstruction using simulated annealing. 

a) Original, Initialization and Reconstructed objects 

b) Evolution of the criterion : J = Ji + AJ2 where Jl — Q{z) and J2 = il{z) 

In Fig. ^ we show a comparison between the results obtained by the proposed 
method and those obtained either by a classical backprojection method or by some 
other methods in the first approach using (||) and (jij) with different regularization 
functionals ^l{f) among those in (||). Also, for the purpose of curiosity we show 
the result of a binary segmented image obtained by thresholding these last images. 

4. Conclusions 

A new method for tomographic image reconstruction of a compact object from its 
limited angle projections is proposed. The basic idea of the proposed method is 
to model the object as a polygonal disc whose vertices coordinates are estimated 
directly from the projections using the Bayesian MAP estimation framework or 
equivalently by optimizing a regularized criterion. 

This criterion is not unimodal. To optimize it two methods are examined: a 
global optimization method based on simulated annealing and a local gradient- 
based method with a good initialization obtained using a moment based method. 
The first one seems to give entire satisfaction and better results. The final destina- 
tion of the proposed method is for non destructive testing (NDT) and evaluation 
(NDE) image reconstruction applications including X-rays, ultrasound or Eddy 
currents [|9|, 
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o:orig, +: init, *: rec. ^^q* ciitere J=J1-kI2 




-60 -40 -20 20 40 60 2 4 6 8 10 12 14 



Figure 4: Reconstruction using a moment-based initialization and a local mini- 
mizer. 

a) Original, Initialization and Reconstructed objects 

b) Evolution of the criterion J = Ji + XJ2 during the iterations. 

References 

1. G. Herman, H. Tuy, H. Langcnbcrg, and P. Sabatier, Basic Methods of Tomography 
and Inverse Problems. Adams Hilger, 1987. 

2. A. Kak and M. Slaney, Principles of Computerized Tomographic Imaging. New York, 
NY: IEEE Press, 1987. 

3. S. Geman and D. McClure, "Statistical methods for tomographic image reconstruc- 
tion," in Proc. of the 46-th Session of the ISI, Bulletin of the ISI, vol. 52, pp. 22-26, 
1987. 

4. G. Demoment, "Image reconstruction and restoration : Overview of common es- 
timation structure and problems," ieeeASSP, vol. ASSP-37, pp. 2024-2036, Dec. 
1989. 

5. S. Brette, J. Idier, and A. Moiiammad-Djafari, "Scale invariant Markov models for 
linear inverse problems," in Proc. of the Section on Bayesian Statistical Sciences, 
(Alicante, Spain), pp. 266-270, American Statistical Association, 1994. 

6. A. Mohammad-Djafari and J. Idier, A scale invariant Bayesian method to solve 
linear inverse problems, pp. 121-134. Maximum entropy and Bayesian methods, 
Santa Barbara, U.S.A.: Kluwer Academic Publ., g. heidbreder ed., 1996. 

7. C. Bouman and K. Sauer, "A generalized Gaussian image model for edge-preserving 
MAP estimation," IEEE Transactions on Image Processing, vol. IP-2, pp. 296-310, 
July 1993. 

8. L. Bedini, I. Gerace, and A. Tonazzini, "A deterministic algorithm for reconstructing 
images with interacting discontinuities," Computer Vision and Graphics and Image 
Processing, vol. 56, pp. 109-123, March 1994. AMD. 

9. M. Nikolova, A. Mohammad-Djafari, and J. Idier, "Inversion of large-support ill- 
conditioned linear operators using a Markov model with a line process," in ICASSP, 
vol. V, (Adelaide, Australia), pp. 357-360, 1994. 

10. M. Nikolova, J. Idier, and A. Mohammad-Djafari, "Inversion of large-support ill- 



A BAYESIAN APPROACH TO SHAPE RECONSTRUCTION 



9 



O'lOirioil Backptojection Gauss Markov 




-60 -40 -20 20 40 60 -60 -40 -20 20 40 60 -60 -40 -20 20 40 60 



b d f 

Figure 5: A comparison with backprojection and some other classical methods 
a) Original, b) Proposed method, 

c) Backprojection, d) Binary threshold of c), 

e) Gaussian Markov Random Field (GMRF) modeling and the MAP estimation 
reconstruction, f) Binary threshold of e). 
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